knitr::opts_chunk$set(
error = FALSE,
message = FALSE,
warning = FALSE
)
library(tidyverse)
library(repurrrsive)
library(sf)
library(wesanderson)
library(leaflet)
library(tmap)
Question 1
wesSelectPals <-
wes_palettes %>%
map_dbl(length) %>%
keep(. > 4) %>%
names()
print(wesSelectPals)
## [1] "BottleRocket1" "BottleRocket2" "Rushmore1" "Rushmore"
## [5] "Royal2" "Zissou1" "Darjeeling1" "Darjeeling2"
## [9] "FantasticFox1" "Moonrise3" "Cavalcanti1" "IsleofDogs1"
## [13] "IsleofDogs2"
Question 2
trade_balance <- read_csv("data/trade_balance.csv")
top_five_countries_names <-
trade_balance %>%
group_by(country) %>%
summarize(balance = mean(balance)) %>%
slice_max(balance, n = 5) %>%
pull(country)
plotWes <- function(palette, df, top_five){
data <-
df %>%
filter(country %in% top_five)
ggplot(data = data, aes(x = year, y = balance, color = country)) +
geom_line() +
scale_color_manual(values = wes_palette(palette)) +
scale_y_continuous(labels = scales::label_comma()) +
labs(
x = "",
y = "Import/Export Balance (Million USD)",
title = "Top 5 Countries by Net Import, 1995-2019",
subtitle = paste("Wes Anderson Palette", palette),
color = ""
) +
theme_light()
}
map(wesSelectPals, ~ plotWes(., trade_balance, top_five_countries_names))













Question 3
counties <- read_sf("data/cb_2018_us_county_20m") %>%
mutate(
fips = paste0(STATEFP, COUNTYFP),
areaCalc = st_area(.) %>% units::set_units("us_survey_mile^2")
)
acs_names <- read_csv("data/ACS2018Counties.csv", n_max = 1) %>% names()
acs <-
read_csv(
"data/ACS2018Counties.csv",
skip = 2,
col_names = acs_names
) %>%
janitor::clean_names()
counties_acs <-
counties %>%
left_join(acs, by = "fips") %>%
st_transform(
crs = "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
) %>%
mutate(
popDens = as.numeric(total_population) / areaCalc
) %>%
select(
fips,
name = qualifying_name,
areaCalc,
total_population,
popDens,
geometry
)
head(counties_acs)
## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -291134.5 ymin: -996650.3 xmax: 1545732 ymax: 343262.2
## CRS: +proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
## # A tibble: 6 × 6
## fips name areaCalc total_population popDens geometry
## <chr> <chr> [us_survey_mile^2] <dbl> [1/us_survey_mile^2] <MULTIPOLYGON [m]>
## 1 37017 Bladen County, North Carolina 880. 33778 38.4 (((1475757 -473449.3, 14…
## 2 37167 Stanly County, North Carolina 405. 61114 151. (((1332365 -454797.5, 13…
## 3 39153 Summit County, Ohio 423. 541810 1280. (((1125940 219233.4, 112…
## 4 42113 Sullivan County, Pennsylvania 453. 6177 13.6 (((1493945 339398.3, 154…
## 5 48459 Upshur County, Texas 590. 40769 69.0 (((75791.7 -860131.9, 75…
## 6 48049 Brown County, Texas 952. 37834 39.7 (((-288114.9 -923209.2, …
Question 4
quantile <- colorQuantile("OrRd", counties_acs$popDens, n = 5)
leaflet() %>%
addProviderTiles(provider = "OpenStreetMap.Mapnik") %>%
addPolygons(
data = counties_acs %>% st_transform(crs = 4326),
weight = 1,
color = ~quantile(popDens),
fillOpacity = 1
) %>%
setView(lng = -120, lat = 52, zoom = 2.7)
Question 5
ca_air_quality <-
read_csv("data/ad_viz_plotval_data.csv") %>%
janitor::clean_names() %>%
mutate(
statefp = "06",
fips = paste0(statefp, county_code)
) %>%
group_by(fips, county) %>%
summarize(maxpm25 = max(daily_mean_pm2_5_concentration))
counties_ca <-
counties_acs %>%
inner_join(ca_air_quality, by = "fips")
counties_ca_centroid <-
counties_ca %>%
st_centroid()
Question 6
counties_ca_sp <- as(counties_ca, "Spatial")
# REDO BY DROPPING N TO REDUCE RESOLUTION
# SET GRID
ca_grid <-
st_make_grid(counties_ca, n = 50) %>%
st_transform(crs = st_crs(counties_ca))
# I had to stop the assignment at this point - the IDW operation ran for over 3
# hours and did not complete.
idw <- gstat::idw(
formula = maxpm25 ~ 1,
counties_ca_sp,
newdata = ca_grid,
idp = 2
)
## [inverse distance weighted interpolation]
idw <-
idw %>%
st_transform(crs = st_crs(ca_grid))
aqRaster <- stars::st_rasterize(idw["var1.pred"])
tmap_mode("plot")
tm_shape(aqRaster) +
tm_raster(title = "Max Daily PM2.5") +
tm_layout(legend.position = c("right", "top"), frame = FALSE) +
tm_shape(counties_ca) +
tm_borders() +
tm_shape(counties_ca_centroid) +
tm_dots()

Question 7
aqRaster_sf <-
aqRaster %>%
st_as_sf() %>%
st_transform(crs = st_crs(ca_grid))
counties_ca_aq <-
counties_ca %>%
st_intersection(aqRaster_sf) %>%
group_by(fips) %>%
mutate(
area = as.numeric(areaCalc),
weight = area / sum(area)
) %>%
summarize(avg_aq = sum(weight * var1.pred))
ggplot() +
geom_sf(data = counties_ca_aq %>% st_transform(crs = 4326), aes(fill = avg_aq)) +
labs(
title = "California Air Quality by County",
subtitle = "Average Daily Concentration of PM2.5 particles",
fill = "Avg. daily PM2.5"
) +
scale_fill_gradient(low = "lightyellow", high = "darkred") +
theme_void()
